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Abstract 

We investigate the evolution of primordial black hole mass spectrum by including 
i . both accretion of radiation and Hawking evaporation within Brans-Dicke cosmology. 

£^ We also consider the effect of evaporation of primordial black holes on the expansion 

dynamics of the universe. The analytic solutions describing the energy density of 
the black holes in equilibrium with radiation are presented. We demonstrate that 
these solutions act as attractors for the system ensuring stability for both linear 
and nonlinear situations. We show, however, that inclusion of accretion of radiation 
delays the onset of this equilibrium in both radiation and matter-dominated eras. 

PACS numbers : 98.80.Cq, 97.60.Lf, 04.70.Dy 

Key words : Brans-Dicke theory, Primordial black hole mass spectrum, accretion of 
radiation. 

1 Introduction 

The Brans-Dicke (BD) theory p], a scalar-tensor theory of gravity, is regarded as a viable 
alternative of Einstein's general theory of relativity (GTR). In the BD theory, gravitational 
constant G is replaced by a time- dependent scalar field which couples to gravity with a 
coupling parameter oj. BD theory goes over to GTR in the limit u — > oo. Interestingly, 
BD type models arise as low energy effective actions of several higher dimensional Kaluza- 
Klein and string theories [2|. BD theory has been used to understand many cosmological 
phenomena such as inflation [3], early and late time behaviour of the universe [I], cosmic 
acceleration and structure formation [5] and the coincidence problem [B] etc. 

The primordial black holes (PBHs) are the black holes which could be formed in the 
early universe in wide mass range through various mechanisms such as inflation [7J |8], 
initial inhomogeneities [9j [10], phase transition and critical phenomena in gravitational 
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collapse [TU [JS], bubble collision [13], or the decay of cosmic loops [14] . The formation 
masses of PBHs are so small that some of them are completely evaporated by the present 
epoch due to Hawking evaporation [15]. PBHs which evaporate in early times could 
account for baryogenesis [TBI [TT] in the universe. On the other hand the longer lived PBHs 
could act as seeds for structure formation [18] and could also form a significant component 
of dark matter [19] . The most relevant issue concerning PBHs is their longivity which 
depends crucially on the effectiveness of various accretion processes. 

In BD theory the possibility of black hole solutions was first proposed by Hawking 
[20] . Using scalar-tensor gravity theories Barrow and Carr [21] studied PBH evaporation 
during various eras. It has been recently observed that in the context of Brans-Dicke 
theory and it's generalised version, inclusion of the effect of accretion of radiation leads 
to the prolongation of PBH lifetimes [221 123] ■ 

From the very beginning of the discovery of Hawking evaporation, it has been appre- 
ciated that the presence of a population of primordial black holes can have significant 
consequences for the evolution of the universe [24] , which further motivates a discussion 
of the expected mass spectrum of primordial black holes. A number of works [25] have 
been done in this regard for studying the evolution of a continuous mass spectrum of black 
holes. In this context, it is assumed by Carr [2S] that black holes exist with a power law 
mass spectrum. He also pointed out that formation of black holes over an extended mass 
range is possible only if the density perturbation spectrum in the universe has a power 
law Harrison-Zel'dovich form. A particular situation where primordial black holes form 
with a narrow range of masses was studied in the Reference [27] . 

A detailed study on the evolution of population of primordial black holes in the very 
early universe was done by Barrow et al. [28]. Where they included the effect of evapora- 
tion of PBHs on the expansion dynamics. We extend this analysis by including accretion 
of radiation by PBHs and taking Brans-Dicke theory as the theory of gravity in our 
present paper. We use a complete set of network equations which follow the evolution of 
a population of primordial black holes in an expanding universe containing either matter 
or radiation. The complete evolution equations consist of the evolution of the black hole 
spectrum, Friedmann equation describing the expansion of the space-time and an energy 
equation governing the evolution of the radiation energy density. 

2 Evolution equations for the spectrum of black holes 
in Brans-Dicke theory 

Considering both accretion of radiation and Hawking evaporation, we obtain the equations 
which will determine the evolution of the spectrum of black holes and the effect of black 
holes and radiation on the evolution of the scale factor. To begin with, let us assume an 
isotropic Friedmann-Robertson- Walker space-time of zero curvature, i.e. 

ds 2 = -dt 2 + a 2 (t)(dx 2 + dy 2 + dz 2 ) (1) 

where a(t) represents the cosmic scale factor. 



2 



Assuming the number density of the initial black hole spectrum by a simple power-law 
form [25], we can write the initial number density of black holes between masses m and 
m + dm as 

N(m, t = 0)dm = Am~ n Q(m - km p i)dm, k > 1 (2) 

where A represents the amplitude of the spectrum and n represents the spectral index. 
We here introduce a minimum cut-off km p i through the theta function for the formation 
mass of the black holes in the spectrum where k is an arbitrary dimensionless constant 
and m v i as the planck mass. The definition of theta (G)-function ensures that no black 
hole is formed with mass less than and equal to the cut-off mass. This initial cut-off 
helps to avoid the divergences at low masses. The above power law is assumed to hold 
for n > 2 for which the total energy density does not diverge at large masses. As we 
take the accretion of radiation by black holes, the black holes take comparatively more 
time for their complete evaporation. During evaporation the mass of the black hole starts 
to decrease so the cut-off is reduced to make the number density non- zero. This cut-off 
vanishes when the initially lightest holes evaporate completely. From this initial spectrum 
we can determine the spectrum at all later times. Then at a given time (t), the number 
density of black holes with cut-off in the spectrum will be 



N(t) = / N(m,t)dm (3) 
Jo 

Hence the energy density of the black holes becomes 

POO 

PBuit) = / N(m,t)mdm (4) 



where the integrals are carried out at constant time. In order to derive analytical results 
we assume all black holes form at the same time. The black hole number spectrum varies 
with time due to dilution by the expansion of the universe and the evaporation of the 
holes into radiation by Hawking process. Because of the Hawking evaporation the rate of 
decrease of the mass rriBH of a single black hole is given by 

m BH = -4irr 2 BH a H TB H (5) 

where tbh = ^Guibh represents the radius of the black hole, an represents the Stefan- 
Boltzmann constant multiplied with number of degree of freedom available for radiation, 
Tbh = sndm BH * s ^ ne Hawking Temperature and G is a time dependent quantity in 
Brans-Dicke theory. Using these expressions equation (jSJ) becomes 

rriBH = 2 (6) 

2567r d G z ~™~ 



m 



BH 

In Brans-Dicke theory, Barrow and Carr [29] have obtained the following solutions for G 
for different eras, as 



G{t) = { y:U (7) 
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where t e represents the era of radiation-matter equality, to represents the present time, 
G represents the present value of G ~ t p i/m p i and n\ is a parameter related to u as 
ni — i+^J- Since solar system observations [30] require that u be large (co > 10 4 ), ri\ is 
very small (m < 0.00007). 

Using the values of G(t) from the equation ([7]) we can write equation ([6]) for radiation- 
dominated era as 

_ __&h 1 (R s 

mBH ~ 256tt 3 G 2 Cf e )^ml H [ ) 

As a black hole accretes in radiation dominated era, its mass increases from its initial 
value (m ) at the time t = to a maximum value m c at the time t = t c [22] , where 

m 



m c = (1 _ 3 ^ ( 9 ) 

and 



2^-f-l a g (t e \3m 
3^ 2567r3Gg U > 



1/2 

(10) 



which demands that f < 2/3. Till the time t c accretion by black hole remains dominant 
over evaporation and after that evaporation becomes dominant. Therefore, considering 
the evaporation alone from t c onwards the mass of the black hole reduces from m c to some 
other value m at any time t. Now we can solve equation ([8]) by integration with proper 
limits and get 

ml M (t) = ml - — — T —— (t -t c ) (11) 
BHK J c 256vr 3 G2(|)2n 1 v V \ ) 

The time at which the black hole is completely evaporated is now found from equation 
(HD as 

2567T 3 ^ 2 /t \ 2 ^ 3 



tevap = ~^ G If) m c + ^ (12) 

We are, however, interested to evaluate the spectrum at time t by taking accretion and 
evaporation into account. Then the black holes having masses between mo and mo+dmo 
evolve to masses m and m + dm in a time t, where the number density at time t remains 
same to the original number density at time t = 0. Now using the Jacobian factor we 
obtain the number density of the black hole spectrum at time t in radiation dominated 
era as 



N(m,t)dm = A[l - -f 



-n+l 



m 



3clh 



256tt 3 G 2 Cf) 2ni 



m-\t - t c ) 



6 



km. 



m 



:i -if) 



3a 



H 



256vr 3 G 2 Cf) 2n i 



f) k-*m-?{t-t c ) 



n + 2 
3 



dm 



(13) 



Equation ([TBI tracks the behaviour of the cut-off which at later time takes the form 



m cut-off(t) 



km 



pi 



./') 



3« 



H 



256tt 3 G 2 Cf) 2 ^ 



l--f) k-*m-?{t-t c ) 



(14) 
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This cut-off found through the theta (G)-function shows that as long as the accretion is 
dominant (t < t c ) the cut-off value increases and beyond t c evaporation becomes dominant 
making the cut-off value monotonically decrease. The cut-off value reaches zero at the 
time 



t 



V256tt : - 



( 



2m 



k 3 m 3 pl (l 



-J 



tr 



(15) 



and vanishes from the picture after that time. The number density of black hole spectrum 
decreases at late times when the cut-off vanishes due to the complete evaporation of 
lightest blackholes. Substituting equation f[T3~j) in equation (jlj) the black hole energy 
density becomes 



PBHif) 
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A(l-- 2 f 



m + I 1 + 



km 



m — 



pi 



:i - if) 



3clh 



256tt 3 Gg(fa)2m 



256tt 3 G^(|) 2 "i 
3 x3 
~ 2 



m-'\t-t c ) 



n + 2 
3 



dm 



(16) 



The amount of energy transferred by the PBHs through evaporation between times t 
and t + dt is given by the mass difference of the black holes at these two times, i.e. 



dE = M(t) - M{t + dt) 



dM 



dt 



dpBH 

dt 



dt 



(17) 



As the cut-off position varies with time we can compute the rate of the energy transfer 
by using equation ( TTB"j) as 



dE 

—r= A 
dt 



n + 2\ 3a H 



3 7 2567r 3 G'o 2 (f a ) 2ni 



3 - f 



-n+l 



m 



-n-2 



m-'\t-t c ) 



-(2±£) 



3clh 



256tt 3 G 2 Cf) 2 ^ 



3 3 

pi 



k m 



dm 



A 



-n-l„ 



m 



m cut-off(t) 

i 3a H 



1 + 



3« 



H 



256vr 3 
3 



P i 



e 



fc 3 mp Z (l 



!/)- 



3ag 
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where 
m cut - ff{t) 



max 



0. 



1-1/) 



3a^ 



256vr 3 Gg(|a)2m 



/) Ar 3 m; 3 (t-t c ) 



(19) 



From equation (|T8|) it is clear that the first term is the only contributor at late times. 
As the cut-off disappears at sufficiently late times the second term vanishes due to the 
presence of O-function in it. 

Now invoking expansion of the universe, as the volume of the universe increases the 
black hole density starts to dilute. Hence the spectral amplitude A will be replaced by 
Aa~ 3 , where a(t) = ^1. From this definition we have ct(0) = 1 at the formation time 
of black holes (t = 0). Now, the equation of motion containing Friedmann equation for 
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describing the expansion of space-time and energy equation governing the evolution of 
the radiation energy density p R in radiation dominated era in Brans-Dicke theory can be 
written as, 

/d\ 2 87T /tn\ ni / \ , N 

(a) = T G "(d <'« + '-> (20 > 

. a dE , , 

p R = -4-p R + — 21 
a at 

In order to solve equations f[2"U]) and f[2~T]) we have to assume sufficiently late times so 
that the lightest black holes are completely evaporated making the cut-off vanish from 

-1 / \ 2ni 



the system. This corresponds to the time t > 
Implementing this condition we can evaluate pbh and ^ as 



(^r)"^o(g) ni ^(l-l/)- 3 + t c 



PBif(t) = Pbho (l - If) ~ n+1 a- 3 (t) (t-t c y , n > 2 (22) 
and 

^ = E ( 1 - If) ~ n+l a-\t) (t-tc)^ 1 , n>2 (23) 



dt 

2 — n 



dt "\ 2 
where Pbh and £7 are constants given by 



A( 3a H 1 Nt 1 /n-2 4\ 
^o--3^ 2567r3G 2 ( ^ 2 nJ 5 v^'3J (24) 



2-n 



Eo -'^(^) (256^ 3 ^(^4) (25) 



where -B^ 21 ^, |j and -B^ 1 ^, |J are the Beta functions. 
From equations ([22"]) and (123]) we find 



dE 1 (n - 2) 



rft 3 (t-t e ) ^ 

Using equations f ]2"2"]) and ([25]) in equations f[2"U]) and f[2"T]) the evolution equations of the 
universe in radiation dominated era can be written in a simpler form as 

« 2 W = yG (^) a 2 (t)[^(t)+p^o(l-|/) a^(t)(t-t c }j ] (27) 

= - 4 ^^W + ^^o(l-^)^ + V 3 (t)(t-t c )^ (28) 

It is clear from equations (|25j) . (|27|) and (|2B| ) that in the limit / -> and t c = 0, 
these equations take the form as obtained in Barrow et al. [28J. It may be noted that 
non-zero / and t c capture the effect of accretion of radiation by the PBHs in radiation- 
dominated era. From equation (|22]) . we find that with increase in both accretion efficiency 
(/) and upper limit of accretion time (t c ) black hole density increases. Similar observations 
apply for the rate of energy transferred by PBHs through evaporation, which is concluded 
from equation (1231) . 
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3 Black holes in a radiation- dominated era 



3.1 An exact equilibrium solution 

We begin by considering a power-law solution of equations (1271) and (128]) for (t > 0) as 

a(t) = a t r , p R (t) = p Ro t s (29) 

By the use of Equation (|29|) . we get from equations (1271) and (128]) that 



-2. 



(30) 



where we have taken (l — yj as constant which strictly holds for t t c . 

The coefficients of power laws a and p# can be found by using equations (T2"9l and 
(13"U|) into the evolution equations f[2"7j) and (T2"g]) as 



87tG V 9 



n\ 2 /t e \ n ± 



2(7 -2n) 
3(n - 2) 



1 - — 



(31) 



PBH Q ot Q ' 



-n) 2 (l -2n) /t \ _ni 



3247rGV 



n-2 



3 n«-i 
2 ; 



1 - - 



t r \ (n+l)/3 



1 + 



2(7- 2n) 
3(n - 2) 



1 - - 

t 



(32) 



where «o is determined in terms of the known Pbh which is given by equation (I24j) . These 
values of coefficients imply that the above power-law solutions exist only for n G (2, 7/2). 
With this range of n, equation (130]) implies r e (1/2, 2/3). Now equation (129]) confirms 
that the expansion power law is intermediate between that of pure radiation-domination 
and nearly matter-domination considering negligible value of ni. Hence the presence of 
black holes results in a different expansion rate compared with that of a conventional 
Friedmann cosmology. 

Now the ratio of energy densities for this solution can be calculated as 



- - 



3 V n - 2 



(33) 



This equation shows that the ratio depends on spectral index n and y. But at sufficient 
late time, y — > and the ratio becomes a fixed constant depending only on the spectral 
index n. So the solution represents an equilibrium between black holes and radiation at 
a late time t decided by time t c beyond which evaporation dominates over accretion of 
radiation. The late time (t > t c ) equilibrium exists because the decrease in radiation 
density due to the expansion of space-time relative to the black holes is compensated by 
the radiation energy evaporated from the black holes. 
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At sufficiently late time the above ratio gives that as n — > 2 from above the black 
hole density dominates and correspondingly r — >■ 2/3 giving nearly matter-dominated 
state as n\ is very small, whereas n — > 7/2 from below radiation density dominates and 
correspondingly r — > 1/2 gives pure radiation-dominated state. This is practically the 
same result as obtained by Barrow et al. [28] though in our case it is true only for t 3> t c 
due to consideration of accretion of radiation. 



3.2 Stability of the equilibrium solution 

It is pointed out by Barrow etal. [28] that the above power-law solution is a particular 
one for its fixed coefficients and fixed powers of t. Such a solution will not be valid for 
arbitrary initial conditions. So they modified the above analysis by using a co-ordinate 
transformation to make the above solution an attractor for a wide class of initial condi- 
tions. Now we undertake similar analysis in BD theory so that both linear and non-linear 
stability of the exact solution can be studied. 

We use equations (1511) and (1321) in equation (1291) to set 



(34) 



Pi?a 4 = pBmV{t)t~ 



(35) 



where £ and rj are two new dynamical variables. Here we replace time t by a new time 
co-ordinate r such that t = e T and we indicate all the derivatives with respect to r by 
primes. Using these equations we can write the evolution equations (12~T1) and (j2"5|) in terms 
of the new variables £ and n as 




— n 



(36) 



14 - An- 

V = -[ n )V + 



n 



1 



3 -j 



-ra+l 



tc 

1 - - 
t 



(37) 



9 / ' V 3 

In these new coordinates, the exact solutions of the above two equations correspond to a 
critical point which is found in the positive quadrant of the phase plane as 



'87^0 



9 



n 



I 



3 \ - 
2. 



1-*= 



— n 



14 - An t 



(38) 



3(n - 2) 
2(7 -2n) 



1 



5/ 

2^ 



-n+l 



1 - - 



'87rG / 9 



2/t 
- // / \t 



2/ 

2^ 



-n+l 



1 - - 



14 - An t 



(39) 
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Above values of £o an d ?7o imply that accretion plays a role in fixing the critical point. 
For increasing / subject to condition / < 2/3, £o an d rj increase. Also accretion delays 
the approach to the critical point. But for a particular value of / and at sufficient late 
time, £ and rjo become constant and act as a true critical point. 

The Sketch of the phase plane for the system (for / = and large t) is shown in figure- 
1. From this figure, we find that all initial conditions tend towards the exact solution at 
late time. 




Figure 1: The phase plane for black holes in equilibrium with radiation in an expanding 
universe. Our analytic solution corresponds to the critical point (£ , Vo) given by equations 
fl38|) and fl39l) and is seen to be an attractor for all initial conditions in the first quadrant 
(the physical region). 



Now we study the linear and non-linear stability of the critical point by following 
analysis of Simmons [31]. First we translate the co-ordinates so as to bring the critical 
point to the origin. Thus our new variables are defined as p = £ — £o and a = r\ — rjo- 

Linearizing with respect to these variables we obtain the linear stability equations as 



P 



1 + n 



P 



7-2n 
9 



2. 



l-f) a 



(40) 



a 



n 



-n+l 



1 - - 



t.\(-"- 1 )/ 3 



3 JV 2 J J V t 
From these, we obtain the auxiliary equation as [31 



P 



14 - 4n\ 
9 J 



(41) 



V 



5 — n 



1 

V H 

y 81 L 



1 ■■■■ ;/)(] I - I//) - 3(// - 2H7- 2//)( 1 - f j 



(42) 
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The roots of this auxiliary equation is given by 

H[- (^W( W - 8> + -)< 14 - 4 -» - 3 <- - 2 >< 7 - *4 - IT}] 

(43) 

But the roots of the auxiliary equation actually determine the stability of the critical point. 
Since the roots are both real, distinct and negative, the critical point is asymptotically 
stable and the same sign (i.e. -ve) implies that it is a node i.e all trajectories enter the 
critical point. 

It may be mentioned here that we had neglected the non-linear terms in setting up 
of the linear set of equation. This stands aposteriori justified as long as one works near 
critical points. Thus the stability shown above for linear system of equations also implies 
stability for full non-linear system. Thus our exact solution is asymptotically stable and 
is an attractor for all physical initial conditions. 



3.3 The approach to the attractor 



Now we estimate the rate of approach to the attractor by using the assumption p BH <C Pr 
in equations (1271) and fl28|) . This assumption actually gives a pure radiation-dominated 
solution as 

^32j[p Ro G ft Q 



a(t) 



Pn(t) 




(44) 



(45) 



The 



Here pr is the radiation density at the formation time of black hole (i.e. t = 0). 
above equation f l4"5]) shows that the radiation density falls off as t~ 2 . 

Again using equation ( 144"|) in equation f )22|) . we can determine the time variation of 
black hole density as 



Pbh = PBH a 3 (t-t c 



>(2-n)/3 



-J 



-n+1 



t 



-(2n+5)/6 



-J 



-n + 1 



t cX (2-n)/3 
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(46) 



Equation (T46]) tells that with increase in the value of accretion efficiency (/) black hole 
energy density increases. But for a constant accretion efficiency (/), black hole density 
is controlled by spectral index (n). For n < 7/2, where actually attractor solution exists, 
black hole density falls at a slower rate than radiation density. Since initial radiation 
density is much larger than initial black hole density, the domination of black hole density 
is possible at a sufficient large time. Closer the value of n to 7/2, more time is required 
for black hole density term to become important in Friedmann equation. 

It is also simple to show that as long as Pbh Pr the energy fed into radiation from 
the holes will never be significant, so there is no possibility that the evaporation can help 
to maintain radiation domination. 
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The approach to the attractor is gauged by the evolution of the quantity 

PBH t ( 7 - 2 n) /6 / 1 _3A--Y 1 _t £ y-)/3 



Pr V VJ V t 

This ratio clearly depends on accretion; with increase in the value of accretion efficiency 
/ and t c , it increases. 

4 Evolution of blackholes in a matter-dominated era 

We now study evolution of black hole spectrum in a matter-dominated universe. For 
this we have to generalise our equations of motion to include matter terms as well as 
radiation terms. Let us introduce an additional parameter /3 to measure the fraction of 
the evaporated energy of the black holes which goes to the radiation and hence (1-/3) is 
the fraction going into matter (i.e. into particles which become rapidly non-relativistic), 
where /3 could in general be a function of the mass distribution of the black holes but here 
it is assumed to be a constant. 

Now the Friedmann equation and energy equations are given by 
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. n 87T /t \ n l 9 T / 3 ,\- n+l o. ,2-n 

a = — G 0[j) OC 2 ^PM + PR + PBH [l--fJ a 6 (t-t r ) ■ 



pr = -i^PR + p(' 1 1 ^)pb«„ (i - §/T" +1 "- 3 e - («) 

a , ^/n — 2\ / 3 „\ - n + l , . . -i-n , , 

Pm = -3-Pm + (1 - 0) [—)Pbho (l - 3/) a (* - ^)~ ( 5 °) 

where /?m represents the matter energy density. Here we have used the assumption that 
we are past the cut-off. Assuming and pbh terms are to be negligible in comparision 
to Pm term in matter-dominated universe we can drop them from the above equations 
and (I50p to get the standard matter dominated solutions to these equations as 



where C is a constant determined from the initial conditions. Now using the above scale 
factor we can obtain an expression for black hole energy density which is given in previous 
section (equation fl2T|) ) as 

/ 3 \ -" + 1 „/ £„\ ^T 1 -4-n + 3ni 

PBH® = PBH (l - ~f) «o 3 (l " j) (53) 
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where Pbh is same as in the previous section and 



J?r^ y/3 ( {2 -m\2 /2-m\ w 2 i-i/3 



<*. = ^ T G„ C j ( + ^j ni - ) (54) 

Now the radiation equation given by equation ( H9|) becomes 

/4m-8\ /n-2\ / 3 _ 3 / t c \^P -7-n+ 3ni 

The general solution of this equation is given by 

n — 2 \ /_ 3 „\ -n+l _ 3 / t„\ "V" -4-w+3n! _ 4w t 



^ = ^13^)^(1-^ V(l"fj f^H-K^" (56) 

where D is an integration constant determined from the initial value of radiation energy 
density. In this way we found the most general solution for the evolution of the black holes 
with matter and radiation in a matter-dominated universe. In this solution for pu, the 
second term gives the conventional evolution of radiation in a matter-dominated universe. 
But at the late times it is the first term which dominates. Now by using equations (153]) 
and (1561) we can calculate the ratio of the energy densities of black holes and radiation as 

P** l ^-n-nA frM (57) 



p R f3\ n — 2 J\ t 

Here again we find for at late times t ^> t c , equilibrium is attained. Thus accretion delays 
the onset of equilibrium. At equilibrium, the above ratio would be a constant depending 
only on the spectral index and ni. Time dependence implies that the black holes and 
radiation tend towards an equilibrium where their energy densities scale the same way 
with time. One may also note that for n > 2 equations ( 152|) . (153]) and ( |56|) imply that black 
hole and radiation densities fall-off faster than matter justifying our earlier assumption. 
So our matter-dominated solution remains self-consistent. 

Proceeding in an analogous analysis as in case of radiation-dominated era one can 
easily establish the attractor property of the equilibrium solutions in matter-dominated 
era aswell. 



5 Discussion and Conclusion 

In this paper we have studied the evolution of a power-law mass spectrum of PBHs by 
including accretion of radiation by PBHs along with Hawking evaporation within the con- 
text of Brans-Dicke theory. Here we also include the effects of the black hole evaporation 
on the expansion dynamics of the universe. We find that in radiation dominated era, 
there exists an equilibrium between energy densities of PBHs and radiation. In order to 
maintain the equilibrium, the PBHs are feeding a substantial amount of energy into the 
radiation via their evaporation. Inclusion of accretion of radiation delays the onset of this 
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equilibrium but once it is set up, accretion plays no further role. The positiveness of the 
ratio between energy densities of PBHs and radiation constrains the spectral index n to lie 
between 2 and 7/2. But we find a stable exact solution in which the scale factor expands 
as a(t) oc t( 8 -™)/ 9 . Thus we see that the scale factor evolves at a rate intermediate between 
that of a pure radiation-dominated (when n — > 7/2) and nearly matter-dominated (when 
n — > 2) universe. Then we study the nature of the stability of the equilibrium solution 
and get that the solution is asymptotically stable critical point of the full system and is 
a node. So for all physical initial conditions, the equilibrium solution act as an attractor. 
We have also shown similar results to hold good in matter-dominated era. 
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